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Abstract — Shannon in his 1949 paper suggested the use of 
derivatives to increase the W*T product of the sampled signal. 
Use of derivatives enables improved reconstruction particularly 
in the case of non-uniformly sampled signals. An FM-AM 
representation for Lagrange/Hermite type interpolation and a 
reconstruction technique are discussed. The representation using 
a product of a polynomial and exponential of a polynomial is 
extensible to two dimensions. 

When the directly available information is inadequate, esti- 
mation of the signal and its derivative based on the correlation 
characteristics of Gaussian filtered noise has been studied. This 
requires computation of incomplete normal integrals. Reduction 
methods for reducing multivariate normal variables include 
multistage partitioning, dynamic path integral and Hermite 
expansion for computing the probability integrals necessary for 
estimating the mean of the signal and its derivative at points 
intermediate between zero or threshold crossings. The signals 
and their derivatives as measured or estimated are utilized to 
reconstruct the signal at a desired sampling rate. 

I. Introduction 

The commonest interpolator is a Lagrange polynomial in- 
terpolator. Widely used Whittacker-Kotelnikov-Shannon [1]- 
[3] interpolator for uniform sampling has a close relation 
to Lagrange interpolation. Shannon in his 1949 paper [3] 
pointed out the possible application of derivatives of a signal 
to increase the WT product. The usefulness of derivatives in 
telemetry was discussed in 1955 by Fogel [4]. The extension to 
non-uniform sampling was developed by Linden and Abram- 
son [5] and Rawn [6]. Interestingly the theoretical framework 
for interpolation using a function and its derivatives was built 
by Hermite more than 130 years ago [7], [8]. A very large 
literature on interpolation and reconstruction now exists [9], 
[10]. Importance of timing accuracy in sampling has long 
been recognized (Papoulis) [11]. This requires greater attention 
when derivatives are used [12]. 

The present work is concerned with a development which 
simplifies the computation involved in incorporating the 
derivatives. The classes of signals considered include natural 
sampling based on threshold crossing and sampling at the 
extrema. Methods for multivariate incomplete integration to 
estimate signal values from correlation characteristic of a 
filtered Gaussian process [13] has been studied. 

Section [XT] discusses the procedure for restoring local sym- 
metry in non-uniform sampling and consequences thereof in 
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simplifying the procedure for incorporating derivative infor- 
mations. Section [HI] is concerned with the framework for 
estimating the statistical mean of the signal and its derivatives 
at a desired time from a knowledge of the correlation structure. 
The techniques of partitioning the correlation matrix or its 
inverse are discussed. Attention is drawn to a path integral 
method in the time domain. This is based on the work of 
Plackett [14], [15]. Hermite expansion [16] for computing 
probability integrals when direct integration proves difficult is 
also considered. Results and discussion are presented in Sec 

ED 

II. Interpolation for non-uniform sampling 

The sine function used in WKS interpolation of uniformly 
sampled signals is symmetric. Chebychev polynomial interpo- 
lation uses non-uniformly spaced zeros but is symmetric about 
the centre. A consequence of non-uniform sampling is that the 
odd derivatives of the polynomial defined by zero locations 
are non-zero. It is useful to locally restore the even symmetry 
about the sampling point. The first derivative of the function 

Go(x) = Y[(l-x/a n )*(l + x/b n ) (1) 

can be removed by multiplying the product by exp(e?l * x) to 
derive 

Gx{x) = Y[(l - x/a n ) * (1 + x/b„) * exp(dl * x) (2) 

where a n and b n give locations of zeros to the right and left 
respectively of the origin and dl = ^(l/a„ — l/b n ). 

More generally the product function Y\(l — x/a n )(\+x/b n ) 
is multiplied by a symmetrizer 

S(x) = exp(dl *x + d3* x 3 /3 + ...) (3) 

where dk = l/ajj — l/&^) for k odd, i.e., to obtain G(x). 
Thus 

G(x) = S(x) * G (x) (4) 

It is to be noted that the even derivatives are necessarily 
non-zero. 

One gets for the case when the first derivative (/'(0)) alone 
is to be incorporated 

f{x) = /(0) * ex P (/'(0)//(0) * x) * G(x) (5) 
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Let f(x) = A(x) * G(x) 



(6) 



where A(x) is the amplitude modulation function and G(x) is 
the switching function or FM term for non-uniform sampling. 
The first few derivatives of A{x) at x — can be found from 
the equations below if dG/dx = and G(x) = 1 at x = as 
desired for any interpolator: 

dA/dx = df/dx (7a) 
d 2 A/dx 2 = d 2 f/dx 2 -3* dA/dx* d 2 G/dx 2 (7b) 
d 3 A/dx 3 = d 3 f/dx 3 - 6 * d 2 A/dx 2 * d 2 G/dx 2 (7c) 

Higher order derivatives of A(x) requires a knowledge of 
lower order derivatives of A(x) and even order derivatives of 
G(x). A formal relation between A(x) and f(x) is derived 
from the expression for the derivative of f{x)/G(x). 

For the case of zero crossing the first derivative f'(x) is 
expressed as f'(x) = A(x) * G(x) and second and higher 
derivatives are derived in the manner indicated. Extremum 
sampling is based on the amplitude and second and higher 
derivatives at points where the first derivative vanishes and 
Eqns. © apply. 

A useful alternative expression when f(0) is not close to 
zero is 



f(x)=f(0)exp{m(x))G(x) 



(8) 



Letting f\(x) = exp(m(x))*G(x), modulation function m{x) 
is derived from the logarithmic derivative of fi(x)/G(x), 
where fi(x) = f(x)/f(0). The above can be stated formally 
as: derivatives of the exponential amplitude modulation m{x) 
are given by the relation 

^rn(x) = ^ (in (m)-^ {H G { x)) (9) 
dx n w dx n \ \f(0) J J dx nK K K " 

The second term in the R.H.S. of Eqn. (O is simply related 
to dn. f(x) given by Eqn. ([8]) is seen to be a product of a 
polynomial and exponential of a polynomial determined by the 
derivatives of the signal. G(x) can be raised to a desired power 
as in Hermite interpolation. In polynomial based generalized 
Hermite interpolation, the order of the polynomial for specified 
zero location is strictly related to the number of derivatives 
desired. This is relaxed in envelop FM description. It is to be 
noted that this operation reduces the contribution from samples 
distant from the point examined, thus reducing as expected 
the number of sample points. Imposition of local symmetry is 
therefore especially useful when derivatives are used. 

A representation of entire functions as a product of a 
polynomial with specified zeros and an exponential of a 
polynomial is useful for nonuniform sampling. This ensures 
that zeros continue to occur at locations desired while the AM 
envelope is determined by the derivatives. 

In the case of a band pass signal with in-phase and 
quadrature components I and Q, one may separately find the 
interpolated values and later combine to form a complex signal 
at a desired frequency. 



In the two dimensional case, one may express 

f(x, y) = /(0, 0) * exp(m(x, y)) * G{x, y) (10) 
In separable form 

G{x, y) = X(x)Y(y) (11) 
Differentials in Eqn. (0 are now replaced by two dimen- 



sional derivatives, i.e., 

gr+s 

: (m(x,y)) = 



dx r dy- 



Qr+s 

dx r dy s 



In 



f(x,y) 



/(o,o), 

& lnx{x) -^ lnY{y) (12) 

A limitation of the exponential representation is the require- 
ment that the signal amplitude is not close to zero. This is 
avoidable by choice of the crossing threshold. 

Symmetrizer defined by Eqn. ([3]) ensures local symmetry of 
the contribution of the signal at z = 0. Approximate symmetry 
for a wider range restricted to narrowband applications can be 
established by introducing a time shift as given by Lomb [17]. 

Taking the simplest case of two point interpolation, one 
finds that a cubic interpolation requires a knowledge of the 
sample value and first derivative at end points as indicated by 

(x - a) 2 {B0 + Bl(x - b)) + (x - bf.{AQ + Bl(x - a)) 

If one uses linear interpolation, four sample points are nec- 
essary. A general result stated in Davis [8] is: the polynomial 



A k 



p(x)=(x-ar^ -^{x-bf + {x-bT^ ^{x-*T d3) 

with A k = Jj^[f(x)/(x-b) n ] and B k = £,[f( x )/( x -a) n ] 
satisfies the condition that the derivatives of p(x) agree with 
the derivatives of f(x) at a and b. For the case of nearly 
sinusoidal signals defined by zeros and specified slopes s\ 
and S2, one may express the function as f(x) = sin(x)(s2X + 
(1 — x)si). Use of second derivative enables one to represent 
functions with two maxima and a minimum as shown in Fig. 
Q] with polynomial only and sin(x) multiplied by exponential 
of a polynomial. 

III. Estimation of signal from correlation 

CHARACTERISTICS 

We restrict our attention in this section to time domain 
signals and the symbols are chosen accordingly. The basic 
assumption of the work of section [Til] is the presence of an 
underlying filter. Linear interpolation over a large number of 
sampling points gives rise to a sine impulse response. When 
the number of sampling points is small, the use of the filter 
response if known enables good recovery. A stand alone two 
point (0, T) interpolation built on the above basis may be 
expressed as 

x{t) = [x(0)m(t) + u{0) mi (t) + w(0)m 2 (i)} G {T) + 
[x(T)m{T-t) + u{T) mi {T-t) + w{T)m 2 {T-i)] G(t) (14) 

where x(0), x(T), u(0), u(T),w(0), w(T) are the values of the 
amplitudes and the first and second derivatives respectively at 
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- Fifth degree polynomial only 

- sin (t) * exponential of a polynomial 



Fig. 1: Waveform for a polynomial only and sin(x) multiplied 
by a polynomial 
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(16) 
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partial correlation is defined as 

"712 - "7i 3 m 23 



■m\2. 3 



V(l-mf 3 )(l -ml 3 ) 



(17) 



and other terms 77713.2 and 77723.1 are obtained by cyclical 
rotation. 

The characteristic function corresponding to Eqn. ( TBI ) is 



t = and t = T and Gq(T) and G(t) are window functions 
which ensure that the individual sample values are not affected 
by the presence of other samples. 777(7;) and mi(t) are derived 
from the filter impulse response. One can include higher order 
derivatives if these are precisely known. 

The development in section [XT] assumes that the values of the 
function and its derivatives are known. The first derivative at a 
sampling point may not be difficult to measure. Higher order 
derivatives even when measured are likely to be contaminated 
with noise. 

In natural sampling where sampling instants are determined 
by crossings of specified threshold the spacing between two 
successive samples may be wider than the Nyquist interval 
even when near equivalence of derivatives to additional sam- 
ples is assumed. 

A large class of signals belongs to sampled values of a 
filtered Gaussian process with a specified correlation function 
and processes derived from the Gaussian. We assume that 
in the interval of interpolation a few sampled values and 
their derivatives are known and a few more are desired to 
be estimated with some measure of reliability as in the case 
of recovery of missing signals [9]. 



A. Partitioning and reduction of order of probability integral 

A fc— variate normal distribution for a vector x with a 
covariance matrix, M with elements 777.^, is expressed as 



<&(l)) = exp(— l/2(m,ij(jJi(jjj)) 



(18) 



W{x\, ....,x k ) 



cxp (-i (x*Ax)) 
(27r) fc / 2 VlM| 



(15) 



where the matrix A with elements a rs is given by A = 
(M) _1 . The elements of A of order three quoted in many 
papers are for ma = 1 



Three basic probability integrals which must be computed 
include the probability integral P = J W(x)dx; the set of 
means given by 777 r = J x k * W(x)dx and product moments 
rij = ave(xiXj) = J XiXjW(x)dx. A measure of the variance 
of the estimated mean may also be required. It is known that 
the product moment ry can be found by differentiating the 
probability integral P with respect to ay. The 777 r s are simply 
related to the probability integral and M. 

This work is concerned with threshold crossings and the ef- 
fect of derivatives and side information at neighboring points. 
The specific problems to be studied concern estimation of (a) 
two point data of amplitude, the first and second derivative and 
(b) four point data of amplitudes and first derivatives. The 
amplitudes are known but the information about derivatives 
may be confined to their correlation behavior and the polari- 
ties. A feature which makes the algebra a bit involved is the 
requirement that the crossing level is ordinarily not zero. Rice 
[13] in his celebrated paper computed the distribution of time 
interval between crossings of a level from a knowledge of joint 
distribution of the variables a;i(r;i),u(= dx±/dt),X2(t2),v(= 
dx2/dt) and evaluating J J uvW(xi,u,X2,v)dudv. An elab- 
oration which uses the same framework is employed. 

The variables are designated as x(ti), x(t 2 ), the first deriva- 
tive u and second derivative w 2 at t 2 , £(£3), the first derivative 
v and second derivative 7/J3 at 7j3,a;(7j4) and x(t) at a point 
intermediate between ^ an d ^3- Their number including the 
amplitude at the point of estimation for the case of two 
derivatives is nine. The nine-variable density function is parti- 
tioned into two sets, one consisting of four amplitude variables 
a;(r;i), ^(7:2), ^(7:3) and xft^) and other consisting of the vari- 
ables including two first derivatives, two second derivatives 
and x(t) at the time of estimation. If the slopes u and v are 
measured, the problem of estimation simplifies considerably as 
one is then required to find conditional distribution of wi, w 2 
and x(t). If polarity alone of u and v are known, one has to 
consider the probability distribution W5 (w, v, Wi, w 2 ,x(t)). 
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1 ) Partitioning for preprocessing: The variables in the 
normal density function are partitioned into two classes: (a)x a , 
those for which the values are known and (b) Xb, those for 
which the values are not known or the polarity alone is known. 

Given the correlation matrix M(M U) Mi 2 ; M 2 i, M22) 
and its inverse A(An, A12; A21, A22) one can rewrite 
W(x) as the product of W(x a ) and VK(xb/x a ) where 

W^(x a ) =exp(-ix*(M i r 1 ).x a )/(2 7 r) m / 2 ^MLj (19) 
and 

W(x h /x a ) = CXp[-i(x b - R a X a )*A 22 (x b - R a x a ). 

(27r)("- m )/ 2 y|A^|] (20) 

where R a = M 2 i(M 1 " 1 1 ) and |A 22 | = |Mn|/|M|. Eqn. © 
shows that integration regime for Xb is modified due to the 
shifts caused by x a . The informations contained in W(x a ) 
influence in two distinct ways. The first is to increase the 
order of the correlation matrix and the second is to introduce 
effective signals S represented by M 2 i (M 1 ^ 1 1 )x a . One can 
rewrite Eqn. (|20T > as 

^(x b /x a ) - C. exp(-l/2.(x b - S) ( A 22 (x b - S)) (21) 

This has an equivalence in chf which is useful when one 
utilizes Hermite expansion. 

A simple useful example of partitioning is computation of 
time interval between upward crossing at x = h followed by 
a downward crossing at threshold k using the method due to 
(Rice) for finding crossing time in the two point case. The 
joint distribution of two space variables with correlation m 
and slopes thereat may be written in partitioned form as 

Wa(xx,x 2 ,m).W(u 1 ,u 2 / (x 1 ,x 2 )) 

where W a ( Xl ,x 2 ) = exp (-A^M^x)) /(2tt.|Mii|) and 
W b (u,S) = exp(-(u- S)'A 22 .(u- S)).|A 22 |/2tt where 
A 22 = (M 22 - M2i(Mii)~ 1 .Mi 2 ) _1 
The signals S3 and S4 are: 



R(t) =ave(uv)(u > 0, v < 0) = 77 + sin lr i)+ V 1 _r i ) 

where r\ is the correlation between u and v. The effect of 
neighboring sampling points can be studied by enlarging the 
sub-matrix Mn and finding equivalent signals controlling the 
timing window. 

The conditional chf corresponding to ( f20b may be written 

as 

<5>{uj b /x a ) = $(w 6 ) exp(juj b .S) (22) 

where uib is the transform variable associated with x b . 

The conditional probability density function W(xb/x a ) 
of dimension five for the two point interpolation problem 
mentioned above may be written as 

W(x b / X a ) = W 5 (ui ,Wi,U2,W2,x(t)/xi,X 2 ,X3,X4,) (23) 

W5 yields the conditional expected value of x(t) on necessary 
integration. 

More generally let L , K , and Af be the number of such 
amplitude, slope and second derivative values that are known 
by measurement and let L equal number of unknown ampli- 
tude variables which constrain the polarity, K equal number of 
slope variables for which only the polarity is specified and M 
that of number of variables for second derivatives for which the 
polarity alone is specified and N be the number of variables 
for which average value has to be found. 

The dimension of the probability integral is then 

k = L + L + K + K + M + M 

The fc— dimensional density function can be reduced by uti- 
lizing the knowledge of Lq + Kq + Mo variables to obtain 
a density function of dimension n = L + K + M. The 
conditional n— dimensional probability density function yields 
the expected amplitude and derivative values by appropriate 
integration. Following the first partitioning to introduce the 
known variables, further partitioning of W(xb) has to be 
carried out. 



«3 



Si 



(m 13 - m 2Z m 12 )ai + (m 23 - m 13 m 12 )a 2 
(l-m? 2 ) 5 

(mi4 — m 2 ±m 12 )a\ + (m 24 — m 14: m 12 )a 2 
(l-mf 2 ) 



The two dimensional integral J uvW(u, v)dudv can be 
found by integrating B(h, k; p). When h and k are both posi- 
tive, one notices that in the region of interest for crossing(u > 
0, v < 0), S3 is negative thus cutting off low positive velocity 
components. Similar remark applies for the negative velocity. 
The end result is that the integration interval is reduced and 
therefore the time window becomes smaller. The probability 
of crossing for positive threshold is also small because of 
dependence on threshold. 

When h = k = 0, one has the well known formula 



B. Statistical mean 

A useful general result is that the statistical means of n- 
variables may be expressed in terms of n probability integrals 
of (n — 1) variables. 

Expressing W^(x) as 



W{y) =C.exp(-iQ(x) 
where Q(x) = ^ and C = 



(24) 
one gets 



dW 

dx r 



(2ir)"/ 2 A /|M| 

{^a rs x)jW (25) 



Integration of equation ( [251 ) yields ^a rs m s , where m s is 
the statistical mean of the s— th variable, as P r . P r is the 
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integral of n — 1 variables with x r replaced by its value at the 
lower limit. 



, 9Q ( 1 , 
P r = C I — — — exp I — —Q ] dx\dx2 ■ ■ ■ dx v 



= C 



dexp(-Q) 
dx r 



dx 



(26a) 



El exp(— hQ r ) , 
(a r ixmi + a r 2xmi + ..) = C \ — — - — /n ;i : = dx (26b) 



(2tt)™/ 2 vW1 

where Q r is obtained by replacing x r in Q by a r and 
integration of RHS is carried out over (n — 1) variables. 

The LHS of ( |26bt is a weighted sum of means. 

The set of simultaneous equations for the vector 
X m [xmi, . . .,xm n ] 



A.X m = (P 1 ,P 2 ,P 3 ,...,P n 



has the solution 



X m = M.(P) 



(27) 



(28) 



where P is the vector (Pi, P2, P3, . . . , P n )- Derivation of Eqn. 
(l27l i assumes that if Q(x) contains a linear term, the variables 
are transformed to remove it. When the effect of the signals 
S as in Eqn. (I2TI 1 is included, Eqn. d28l > becomes 



X m = M.(P) + S 



(29) 



One can associate a signal flow diagram with the above 
equation. It is necessary to note that M is the correlation 
matrix of the partitioned variables. 

The first moment for a correlated pair of variables provides 
the simplest example and is given for the threshold at zero by 



ave(xi) = ave(x2) 
If threshold are at h and k, 



1 



2V2vr 



(1+P) 



, \ Pi + PP2 , P 2 + pPl , ... 

ave(xi) = — / _ and avea;2 — j= (30a) 



2V27T 



2V2tt 



where P x = e -' l "/ 2 erfc ( |= Hp = | and (30b) 



P 2 = e~ k /2 erfc 




(30c) 



It is well known that a two-dimensional probability integral 
can be reduced to a single integral. 



It is convenient to use Owen's [18] procedure for computing 
2D normal integrals using T— functions, where T is defined 

as 

1 /"' „vi,L_, 



T(h,a) = 



1 [ a exv[-\h 2 {l + x 2 )] 



2tt Jo 



1 + x 1 



(32) 



The result is 



B(h.l,: l>] = T[h, l)+T(k,^ 



T h, k - ph _) T[k. H - Pk 



G(h)G{k) (33) 



where G is one dimensional Gaussian integral. 
For the case of three variables, 



ave(xi) 




ave(a; 2 ) 




_ave(a;3)_ 





mn 
m 3 i 



m 12 
m 2 2 
m 3 2 



mi 3 
m 33 



B(ai,x 2 ,x 3 ) 
B(a 2 ,xi,x 3 ) 
B(a 3 ,xi,x 2 ) 



(34) 



where B(ai,Xj,Xk) signifies a two dimensional integral with 
Xj replaced by a,j in Q(x). The result can be expressed in 



terms of partial correlations. When a\ 
gets 



a 2 
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0, one 



ave 0i) = ^371 (| + sin 1 (m2 3 .i) 



+ sin 1 ( m i2. 3 )j + m 13 



sin [m 13 .2. 



(35) 



The first is a self term and the other two are mutual terms 
representing contributions from other two equivalent sources. 
For the case of four sources, four three dimensional inte- 
grals P(ai,x 2 ,x 3 ,Xi), P(x\,a 2 ,x 3 ,Xi), P(x 1 ,x 2 ,a 3 ,x 4l ), 
P{x\,x 2 ,x 3 ,ai) multiplied by the coupling terms generate 
the means. The integrals in the five dimensional case may be 
written as P(a,i,Xj,Xk,xi,x m ) which arise from integrating 
C. exp(— l/2.Q(oj, Xj, Xk, £/, x m ) with Xi replaced by a, in 
Q(x). 

C. Relating a probability integral to a set of integrals of lower 
dimension 

Reduction of probability integral and integral for finding 
mean to a set of integrals of lower order is considered in 
this section. The probability integral, P n {a\, a\, a 3 ..a n ) is 
computed from 



P n (a\,a 2 ,a 3 ..a n ) 



W{x\, x 2 , x 3 ..)dx\dx2dx 3 ..dx r , 



where a\ to a n are the lower limits while the other limit is 
infinity. 

The general method due to Plackett for reduction of order 
may be stated briefly as follows. 



W(x,y,p)dxdy 



exp 



lh Jk 
h 2 - 2h.k.p + k 2 



2(1 -p 2 ) 



.dp/(2wV(l - P 2 ) (31) 



dP 
d(ri2 



W{a 1 ,a 2 ,x 3 ,...x n .)dx 3 ...dx n (36) 



This can be derived directly or from the partial differential 
equation (Plackett 1954) 
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dW 



d 2 W 



d(rij) dxidxj 



(37) 



Use of 071) yields 

dW 

— r = W 2 (h,k,r 12 ) * W n - 2 (x d ,A 2 2,s 1 ,S2) (38) 

o(ri 2 ) 

A method of considerable power is based on conditional 
probabilities (Steck [19]) which works from low order to a 
higher order. A three dimensional integral can be expressed as 
a sum of three two dimensional integrals and four-dimensional 
probability integrals as sum of four three dimensional inte- 
grals. It is well known that computation of such integrals 
for n > 3 poses difficulties. Several methods are available 
for reduction of the probability integrals principally through 
partitioning for the variables into smaller groups in the manner 
stated earlier. 

Path integral: Plackett devised a line integral method for 
finding the probability integral for a point P defined by the 
correlation matrix M once the value at another point K is 
known. Symbolically for i ^ j 



rriij = (1 - t)mij(K) + tm tj (P) 



(39) 



A modification proposed by Pawula is to multiply off- 
diagonal terms by t in the correlation matrix. 

Dynamical path integral: In physical problems concerning 
filtered Gaussian noise where the correlation function and its 
derivatives can be derived, a direct approach is to form the 
time derivative and use the relation 



dW 



uvv \ -> 

~dT = ^ 



di 



dW 



dt dm; 



(40a) 



to derive 



dP(n) s—^ drriij 



*P(n-2) 



dt ^ dt 
Partitioning variables as in d20b . this can be written as 



(40b) 



dP 
~dt 



^rriij J W(ai,a 2 ,A kh x k ,x l )dx k dxi (41) 



where the matrix Aki corresponds to A22 in Eqn. 
For n — 3 one gets 



dP f 

— =m'i2 / W(a 1 ,a 2 ,x 3 )dx 3 + 

m'13 J W(ai,a 3 ,x 2 )dx 2 + rri 23 J W{a 2 , a 3 , x{)dxi (42) 
The (1,2) component of d42l may be written as 

dP U = rn 12 
dt (2tt)3/ Vl-m? 2 " 

/ 1 al-2m 12 aia2+al \ 
CX V{ 2 1=^ ) 

^Jexp(- ^- s f a ™ )dx 3 

where S = (( mi3 ~ m23 "' 12 ) ai +(™ 23 ^ 7 " 13mi2 )) a2 
1 ~ (i-™? 2 ) ' 



(43) 



The integral can be expressed as an error function with 
argument S\. 

If ai = a 2 — a 3 = 0, one has a simple result 




m i3 



m 23 



'13 



(44) 



'23 . 



This can be readily integrated to derive the well known 
expression P 3 = sin ' 1(mi2)+sin ' 1( ^ 3)+sin ' 1(m23)+7r/2 . The 
advantage of the approach lies in providing intermediate values 
of the integrals from small or no correlation to the present 
correlation values. This may be considered to be a backward 
evolution starting with large time separation when correlations 
are negligible. A use of the dynamical path integral is to find 
t^— by computing derivatives of P and a rs and then dividing 
the results. When the time derivative is not known, Plackett's 
method or Pawula's version may be used. A variation is to 
use time values where some correlations are negligible and 
therefore integration is easy to carry out. 



D. Chf based computation - Hermite expansion 

It is well known that multivariable expansion of the char- 
acteristic function and its inversion is useful for evaluating 
probability integrals and statistical means. The starting step 
in this case too is partitioning and including the effect of 
known signals on the uncertain ones. The transform space 
uj is segmented into uj a and to h and transform $ into 
$a(w Q ),$ & (w b ) and $ ab (uj a ,uj b ), i.e., 

Chf(w) = $a(Wa).$b(w b ).$ Qf) (w Q ,W 6 ) (45) 

Integration over u a gives 



$(o>b/x a ) = / $(w) exp(juj a .x a )duja 



(46) 



This is equivalent to Eqn d22l . Hermite expansion is then 
carried out over to h . We consider the specific problem of 
computation of four variable integral. Let a correlation matrix 
have non-diagonal elements defined as m,\ 2 , m\ 3 , mu, m 23 , 

The expansion of the chf contains typical terms 

p q r s t u 
TO 12- TO 13- m 14- m 23- m 24- m 34 . , ,mi, ,m 2 , ,m 3 , ,m 4 
* UJ-, LO^f .Ldo .CO A 



(p!)( g !)(r!)( fl !)(*!)(«!) 

where mi = p + q + r; m 2 = p + s + t; m 3 — q + s + 
u; m 4 = r + t + u. Noting that multiplication by a power of w 
is equivalent to differentiation of that order of the individual 
Gaussian variable, one can find the probability integrals and 
the means. 

Application of Hermite expansion is limited to three vari- 
ables with three off-diagonal terms. Direct four variate expan- 
sion is expensive. The number of terms for the triple product 
is as large as 216. An alternative using bivariate integrals and 
their derivatives is useful. 

chf can be partitioned as 
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where $ a and <£>b are the chfs of the variables xa and xb and 
$ a ;> represents the mutual correlation terms. A four variable 
chf may be written as a product of chfs <l> a and <£>b of two 
pairs of correlated variables and a term $ a h representing 
mutual correlation. The density function is then a product of 
the density functions W a {x\, x 2 , fnyi) and Wb ( x 3> x ii WI34) 
operated upon by terms resulting from expansion of <J> a {,. 4> a {, 
may be written as 



$ Qb = exp(-(mi3Wl£J2+mi4WlW4 + m23W2W3 + m24W2W4)) 

(47) 

The expansion differs from usual single variable expansions 
as it groups a pair of variables and simplifies computation of 
probability integrals. As an example the probability integral 
and the mean for one sided variables for four variates can 
be expressed as power series in mi 2 , TO13, mi4, 77123, TO 34 and 
term by term integration may then be carried out. One may 
also divide the variables into two groups and the mutual 
correlation chf <£> a (, may be expanded as 



p(x) = 



r((m + n)/2) 



m- 13 .TO^ 4 .TO 23 . TO2 4 

(p!)(«!)(r!)(«!) 



(48) 



where 7771 = p + q, 7722 = r + s, 7773 = p + r and 7774 = q + s. 
The probability density function can now be expressed as: 



F a (Dl m \D2 m *)*W 2 (x 1 ,x 2 ,m 12 )* 
F b (D3 m3 ,D4 m *) * W 2 (x 3 , Xi , m 34 ) 



(49) 



where Dl = d/dx l7 D2 = d/dx 2 ,D3 = d/dx 3 ,D4 = 
d/dx^. F a and Fj, represent sums resulting from inversion 
of (l48l i and W 2 are two dimensional integrals. 

E. Comment 

The results of Sec.|III]can be extended for estimating signals 
for some non-Gaussian processes describable as Gaussian 
process with random parameters. A special class is sub- 
Gaussian symmetric alpha-stable process. A stable random 
vector may be expressed as X — \/~{A) * G where G is a 
Gaussian process. If A has a Laplace transform of the form 
exp(— (s)a/2) where s is the transform variable, the chf of 
the sub-Gaussian process is given by 

chfa(w) = exp(-(Qg(w) a /2) 

where Qg is the exponent of the chf of a Gaussian process. 

A non-Gaussian process generated with random scaling of a 
Gaussian process is described in Grigoriu [20]. When the scale 
parameter of variance has an inverted gamma distribution, the 
unconditional density function on averaging over the scale 
parameter becomes 



r(n/2)(77?T/3 2 ) m /2 1 



-(m+n)/2 



1- 



nj3'- 



^M^x 



For n=l,this becomes a multidimensional Cauchy distribution. 

An approach applicable for filtered Poisson process is a sum 
of Gaussians with distinct covariances.lt is known that some 
combinations of amplitude distribution of impulses and the 
filter characteristics result in nearly Gaussian distribution. The 
conditional distribution for a given number of impulses in an 
interval within the filter time window therefore generates a 
Gaussian sum. 

IV. Results and discussion 

It is known that there is an equivalent Nyquist rate for the 
case of non-uniform sampling. The condition for reconstruc- 
tion is that the sampling interval rate lies within a specified 
range of the Nyquist interval. This restriction is relaxed when 
derivatives are available. The Hermite interpolation recon- 
struction procedure requires that the order of the polynomial 
be doubled for a first derivative compared to the order for 
signal amplitude only. A legitimate procedure is to find the 
inter-sampling interval and compute the necessary number of 
derivatives. If this alternative is not implementable, one has the 
option of estimating values of the signal using characteristics 
of signal. 

The signals are therefore a mixture of deterministic and 
partially known components. The number of components be- 
longing to the second category is restricted by the requirement 
of evaluation of probability integrals. As noted in Sec. [ill] 
it is not difficult to compute the statistical mean of five 
variables based on computation of four probability integrals. 
The two point problems with known terminal amplitudes are 
easily solvable. These include finding (a) mean of slopes at 
two points and amplitude at an intermediate point; (b) two 
second derivatives and an amplitude if slopes are known, 
(c) two first and second derivatives each at terminal points 
and an amplitude. For the four point case, the amplitudes 
at four points are known. If the slopes are also known, one 
can estimate four second derivatives and an amplitude at an 
intermediate point. 

The procedure can be stated as follows. The correlation 
matrix and the density function and chf are first partitioned 
for reduction to partially known variables and those to be esti- 
mated. Probability integrals necessary for finding the statistical 
means are then computed. This step requires partitioning and 
reduction to simple one dimensional integrals. 

The known values of the samples and their derivatives 
are then combined with those estimated for an FM-AM 
representation. This is then converted to uniformly sampled 
representation. 

The filter chosen for computation has a correlation char- 
acteristics given by r{t) — exp(—at 2 ).sinc(t). This is used 
to form the correlation matrix of the amplitudes, and the 
first and second derivatives in the manner discussed by 
Rice(BSTJ,1945). The effective signal value is then computed 
for a specified set of sampling points adjoining the region of 
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estimation. The parameter a permits one to cover the behavior 
from monotonic to pure sine filter. In natural sampling the 
sampling event occurs when a threshold is crossed. The time 
interval between two events depends on a combination of filter 
characteristics and local energy. 

Figure |2] shows r(t),dr/dt and d 2 r/dt 2 for a particular 
choice of a — 0.25. The zero of r(t) occurs at that of sinc(t). 
The event pair of an upward crossing followed by a downward 
crossing is related to the correlation between the slopes at time 
interval ti2- The polarity and amplitude of second derivative 
indicate the time window permitted for such events to take 
place. The value of the threshold has an immediate control on 
the actual time window. 

If the thresholds h and k lie on the same side of zero, R(t) 
is compressed as shown by Fig. [3] which assumes h = k and 
local energy parameter is unity. The control of h and k occurs 
through the medium of effective signal mentioned in Sec. Hill 
The monotonic relationship between R(t) and r(t) remains 
unchanged though. 

It is important for computations based on dynamic time 
integral to find how the value of the determinants changes 
with time scale. Once the correlation matrix is partitioned, 
the correlation behavior becomes conditional. The partial 
correlations for a four variate case are plotted in Fig. @] If 
the partial correlations have low values, the computational 
complexity can be reduced. Figure[5]shows the approximations 
to fourth order and third order probability integrals. Figure [6] 
shows the individual path integral terms of the fourth order 
probability integral and their sum is shown in Fig. [7] A final 
time integral upto the time scale of interest yields the desired 
result. 

A knowledge of the correlation characteristics and local 
energy is the basis of signal estimation discussed in sec 3; 
this does not appear to be simple when the time window is 
very short. A straight forward way is to use the data about 
the amplitude and the derivatives near the estimation region 
to interpolate and resample at uniform high rate and then find 
local correlation features. Teager [21] energy operator which 
uses the signal and its first two derivatives offers another 
means. The operator is expressed as 

TK = (dx/dt) 2 - x{t)d 2 x/dt 2 (50a) 
In discrete form it is written as 

TK = x(n) 2 - x(n - 1) * x(n + 1) (50b) 

where x(n) is the nth sample. This can be interpreted as 
r(0) * (1 — m(t)) where m(t) is normalized correlation for 
small time t between x(n — 1) and x(n + 1). A sum of the 
operator values over a large number of samples may be used 
to indicate the second moment of the spectral density. 

The number of crossings at different levels and the time 
interval between crossings can be utilized for deriving infor- 
mations about the energy and correlation. It is known that 
the number of zero-crossings AO is related to d 2 r(t)/dt 2 at 
t = and equivalently to the second moment of the spectrum. 
The crossing rate at a height h is obtained by multiplying 
AO by exp(— h 2 /2a 2 ). A comparison of the crossing rates at 



different heights will therefore give a measure of a. The timing 
distribution at a zero crossing gives a measure of d 2 r(t)/dt 2 
and time between consecutive crossings at h and k furnishes 
a measure of effective signals and therefore dr/dt. 

Statistical averages of amplitude or derivative contain a 
parameter representing corresponding correlation. When the 
agreement between the actual measurements and the estimates 
is satisfactory, one can extract the correlation parameter for the 
time separations involved. 

Concluding remarks; Symmerization of Lagrangian inter- 
polation function is seen to be useful in incorporation of 
derivatives and leads directly to an envelope-FM representa- 
tion. When signal derivatives are not precisely known, the 
values of signal and its derivative at a point intermediate 
between sampling points can be found from a knowledge 
of the correlation characteristic. Two stage partitioning to 
reduce dimensionality of the probability integrals involved, 
once to make use of knowledge of signals and the derivatives 
and then to compute the conditional probability integrals has 
been shown to be an essential tool. Dynamical path integral 
method commends itself when values for a continuously scaled 
time instants are desired. A technique for finding means of 
several variables simultaneously using a set of lower order 
integrals is shown to be useful for statistical estimation of 
Gaussian signals for interpolation. Hermite expansion based 
computation of probability integrals appears to be the simplest; 
the labor involved can be reduced by partitioning. 

The work in this paper made an effort to find a way 
to incorporate partially known signals with various degrees 
of uncertainty in the interpolation process for non-uniform 
sampling. The key issue of finding random signal parameters 
requires further study and empirical confirmation. 
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